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Abstract 

Gravitational instability is one of considerable mechanisms to explain the forma- 
tion of giant planets. We study the gravitational stability for the protoplanetary 
disks around a protostar. The temperature and Toomre's Q-value are calculated by 
assuming local equilibrium between viscous heating and radiative cooling (local ther- 
mal equilibrium). We assume constant a viscosity and use a cooling function with 
realistic opacity. Then, we derive the critical surface density S c that is necessary for 
a disk to become gravitationally unstable as a function of r. This critical surface 
density S c is strongly affected by the temperature dependence of the opacity. At the 
radius r c ~ 20AU, where ices form, the value of S c changes discontinuously by one 
order of magnitude. This S c is determined only by local thermal process and criterion 
of gravitational instability. By comparing a given surface density profile to S c , one 
can discuss the gravitational instability of protoplanetary disks. As an example, we 
discuss the gravitational instability of two semi-analytic models for protoplanetary 
disks. One is the steady state accretion disk, which is realized after the viscous evo- 
lution. The other is the disk that has the same angular momentum distribution with 
its parent cloud core, which corresponds to the disk that has just formed. As a result, 
it is found that the disks tend to become gravitationally unstable for r >r c because 
ices enable the disks to become low temperature. In the region closer to the protostar 
than r c , it is difficult for a typical protoplanetary disk to fragment because of the 
high temperature and the large Coriolis force. From this result, we conclude that the 
fragmentation near the central star is possible but difficult. 

Key words: accretion, accretion disks — instabilities — (stars:) planetary sys- 
tems: protoplanetary disk — (stars:) planetary systems: formation 
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1. Introduction 



Two major models are currently considered for the formation of gas giant planets. One 
is the core accretion model (CA model; Goldreich & Ward 1973; Hayashi 1981; Pollack et al. 
1996), and the other is the gravitational instability model (GI model; Cameron 1978; Durisen 
et al. 2007). In CA model, massive solid cores are made from dust first, and then gas giant 
planets are formed by the gas accretion onto these cores. In GI model, a disk around a protostar 
fragments into pieces by gravitational instability and these pieces become gas giant planets. 
Recently, extrasolar planets are discovered by direct imaging (Marois et al. 2008; Kalas et al. 
2008). These planets are heavier than Jupiter and farther from the central star than Neptune. 
Unless some additional ideas are considered, it is difficult for CA model to explain these planets 
because it is thought that the gas disappears from the disk before the formation of the massive 
solid core (Dodson-Rovinson et al. 2009). On the other hand, GI model has the possibility to 
make gas giant planets far from the central star if the mass of a disk is greater than that of the 
protostar. Thus, GI model attracts attention to explain their formation. 

In order to understand planet formation by GI model, it is desirable to know the physical 
condition for disk fragmentation and the position where the fragmentation occurs. At present, 
two criteria are frequently used to discuss whether a protoplanetary disk is likely to fragment. 
The first is Toomre's stability criterion (Toomre 1964) 



with gravitational constant G, epicyclic frequency n cp , sound speed c s , surface density S, and 
stability parameter Q. Toomre (1964) showed that the infinitesimally thin disk is stable if 
the stability criterion (1) is satisfied. It is necessary to violate this criterion in order for 
fragmentation to occur. The other criterion is Gammie's cooling criterion (Gammie 2001) 



where Q is the angular velocity, t coo i = u(du/dt)~ is the cooling time, u is the specific internal 
energy, du/dt is the total cooling rate, /3 is cooling parameter, and /3 C is the critical cooling 
parameter. Gammie (2001) suggested that rapid cooling is necessary for fragmentation, in 
addition to violating the Toomre criterion (1). Using shearing sheet simulations with constant 
cooling rate and without heating except for artificial viscosity, he showed that a self-gravitating 
disk can fragment only if /3 < 3. If not, disk cannot fragment because the disk can be stabilized 
by internal heating due to the turbulence arising from gravitational instability. This state is 
called "gravitoturbulence" . It is said that the heating owing to the gravitoturbulent dissipation 
automatically controls the Q value close to unity. The gravitoturbulence has the potential to 
stabilize the system where Toomre criterion is violated, while it does not arise in the disk that 
satisfies Toomre criterion. 




(2) 
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The P c obtained by Gammie (2001) is for the specific case (local, 2D, etc.). Many three- 
dimensional simulations were carried out in order to obtain the general value of (3 C (e.g., Rice 
et al. 2005; Clarke et al. 2007; Cossins et al. 2010). However, each simulation suggested the 
different C , and general f3 c has not been clarified yet. Furthermore, it is not obvious whether 
or not P c exists as some value. Meru & Bate (2011a) calculated the disks with different initial 
parameters such as the disk size, the surface density profile, the mass of the disk, and the mass 
of the central star. From their calculations, it is found that (3 C is different among the disks 
with a different set of initial parameters. This result shows that it is not determined only by /3 
whether the disk can fragment or not. Moreover, Meru & Bate (2011b) calculated the /3 C with 
various numerical resolutions. They showed that the disk was able to fragment with larger /3 if 
the resolution is better, and they had no indication of the convergence of (3. From these results, 
it is not clear that Gammie criterion is always applicable to discuss the fragmentation of a disk. 

Some analytical works for disk fragmentation have been done. Rafikov (2005) argued 
that a protoplanetary disk is unlikely to fragment near the protostar (r < 100AU) by using two 
criteria described above. However, as noted above, Gammie criterion is uncertain, and thus 
this argument includes the same uncertainty. On the other hand, Toomre criterion is probably 
assured in the sense that the disk violating Toomre criterion is gravitationally unstable. Other 
previous analytical works for disk fragmentation (e.g., Clarke 2009; Kratter & Murray-Clay 
2011) assumed gravitoturbulent disks with Q — 1. These studies do not have the ability to 
discuss the position where Toomre criterion is satisfied. Since Toomre criterion provides ro- 
bust necessary condition for fragmentation, the studies based solely on Toomre criterion (1) is 
favorable to know the possibility of fragmentation. 

In this study, we investigate the gravitational stability in protoplanetary disks based 
solely on Toomre criterion without the uncertainty associated with Gammie criterion. We do 
not aim to consider the time evolution of particular disks but derive an instability condition that 
is available for various disks. Since recent observation shows the diversity of protoplanetary 
disks and planetary systems, we think that it is important to derive such a widely useful 
formulation. In order to calculate Toomre's Q value with given £ and K cp , it is necessary to know 
the temperature. Thus, it is significant to consider realistic thermal processes. We construct 
an analytical model for temperature calculation with thermal process in section 2. In section 3, 
we derive the critical surface density S c that is necessary for a disk to become gravitationally 
unstable. This S c is available to discuss the possibility and position of gravitational instability, 
regardless of surface density profile or formation process of the disks. In section 4, we introduce 
two semi-analytic models for protoplanetary disks, and discuss the possibility and location of 
the gravitational instability by comparing the surface densities of these models to S c . In section 
5, our results are discussed and compared to previous studies. Finally, we summarize this study 
in section 6. 
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2. Temperature 



In this section, an analytic model for the temperature calculation is introduced. The 
local thermal equilibrium between viscous heating and radiative cooling are assumed. We use 
Keplerian rotation law, ideal gas, vertical hydrostatic equilibrium, and a geometrically-thin 
disk. We assume that radiation escapes only in the vertical direction. Effective viscosity with 
constant a introduced by Shakura & Sunyaev (1973) is adopted. Since these assumptions are 
similar to the standard accretion disk (Pringle 1981), we have 



n = v/GM s /r3, (3) 

Cs = / jk B T/m, (4) 

H = cjn, (5) 

p = £/2# = Sfi/2c s , (6) 

r = -ac s Hm 2 , (7) 



32ctT 4 

A = { 3kS 

8aT 4 K £ 



(t>1) 

(r<r 



and 

r = A, (9) 

where 7 is specific heat, fn is mean molecular weight, p is density, H is scale height, M s is the 
mass of protostar, a is Stefan-Boltzmann constant, k-Q is the Boltzmann constant, V is viscous 
heating rate, k is opacity, A is cooling rate per unit area, and r = Sk/2 is optical depth of 
vertical direction. We approximate cooling rate by connecting optically thin and thick limit 
continuously. Here, we use the Rosseland mean opacity n(p,T) that is approximated by using 
power-law as 

K = K p a T\ (10) 

where the values of k,q, a, b are summarized in table 1 (Bell & Lin 1994; Cossins et al. 2010). 
This table shows the characteristic temperatures at which the main component of opacity 
changes. For example, T\ = 166. 8K is the temperature where the main component of opacity 
changes between ices and sublimation of ices, and T 2 = 202. 6K is between sublimation of ices 
and dust grains. Other characteristic temperatures depend on the density. 

In order to calculate the temperature as a function of r and S, seven variables 
Q, c s , H, p, k, T, A can be eliminated by using the eight equations from (3) to (10). Thus, 
we can obtain the equilibrium temperature analytically for a given r and S. Additionally, 
we introduce the minimum temperature T min in order to mimic the effect of ambient heating 
sources in a molecular cloud. If the above equilibrium temperature becomes lower than T min , 
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Table 1. Bcll&Lin opacity 



Temperature Temperature 



v^jjaeiuy icjiiiiie 


h.o^cm /gj 


Ci 


h 

u 




i aiige _l(jiivi 


Ices 


Z X 1U 


u 


z 


u 


100. o 


Sublimation of ices 


2 x 10 16 





-7 


166.8 


202.6 


Metal dust 


1 x 10" 1 





1/2 


202.6 


2286. 7p 2 / 49 


Sublimation of metal dust 


2 x 10 81 


i 


-24 


2286. 7p 2/49 


2029. 7p 1/81 


Molecules 


1 x 1(T 8 


2/3 


3 


2029. 7/9 1 / 81 


lOOOO.Op 1 / 21 


Hydrogen scattering 


1 x 1(T 36 


1/3 


10 


lOOOO.Op 1 / 21 


31195. 2p 4 / 75 



we simply use T min instead of the equilibrium temperature. As a result, the temperature can 
be represented as 

/ 97 \ V(3-fe+a/2) 

^ Ko A5V- 3 ( a+1 )/ 2 S a+2 ] (r > V 
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From equation (11), we can calculate temperature as a function of £ and r for a given a, 
M s , and T min . In this study, we set a = 0.01, which is the typical value of a caused by MRI 
(magnet rotational instability) turbulence, and assume T min = 10K, which is the typical value 
of a molecular cloud. 

Figure 1 shows the temperature in the r — S plane for the case with M s = 1M Q , a = 
0.01 and T m ; n = 10K as the fiducial case in this study. The solid lines show the contours of 
temperature, and the dotted line shows the r = 1 line. The characteristic temperature 7\ 
defined above is labeled as 166K, and T2 is labeled as 202K. The line for 166K can be regarded 
as snow line because ices become the main component of the opacity in the area below this 
line. Since other characteristic temperatures depend on the density, they do not coincide with 
contour lines of iso-temperature. At these characteristic temperatures, the property of opacity 
drastically changes as shown in equation (10) and table 1. This property of opacity has much 
influence on the gravitational stability of protoplanetary disks. 

From figure 1, it is seen that the temperature is high for large £ and small r. Above 
dotted line, the temperature is high for large surface density because a large amount of matter 
prevents radiation from escaping the system. In a small radius, the dynamical time scale is 
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Fig. 1. Contours of the temperature for the case with M s = 1M Q and a = 0.01 and 
T m i n = 10K. The solid lines show contours of temperature. The dotted line represents r = 1. 



short, and thus the heating rate is large. This is why the temperature is high for large £ and 
small r. The property of cooling rate drastically changes whether or not the condition r > 1 
is satisfied as described in equation (8). Below the dotted line, which means that the disk is 
optically thin, it is seen that the temperature is independent of surface density. This is because 
all radiation emitted from disk can escape in the optically thin case. At the right-bottom area 
of the line with 10K, the temperature becomes the minimum temperature T min due to the small 
viscous heating rate and large cooling rate. We henceforth name this area "isothermal region" . 
Note that the dotted line exists only at low temperature area, such that 10K < T < 20K. This 
is due to the high cooling efficiency for the disk with r ~ 1 . 

3. Critical Surface Density 

Now we can calculate Toomre's Q(r,T,) by using equations (1) and (11). The analytical 
expression for Q value is summarized in appendix 1. In figure 2, the solid line represents 
the critical surface density S c that is necessary to satisfy Q — 1 for a given r. As references, 
the dotted and dashed lines represent the surface densities that satisfy Q = 0.5 and Q — 2, 
respectively. A disk is expected to become gravitationally unstable if the disk has a larger 
surface density than S c . This critical surface density S c is determined only by local thermal 
process and criterion for gravitational instability. Thus, one can discuss the possibility and 
position of gravitational instability regardless of disk formation process by comparing a given 
surface density profile to E c . We show examples of this discussion in the next section. However, 
one should keep in mind some cautions, for example, 1) gravitational instability is not the same 
as fragmentation and 2) Q value is not universal for the case with non-uniform medium or with 
non-axisymmetric situations. 

In figure 2, S c changes discontinuously by one order of magnitude at r = 24 AU. We 
henceforth represent this critical radius as r c . By comparing figure 1 and figure 2, it is confirmed 
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Fig. 2. The critical surface density for the fiducial case. The solid line represents the critical surface 
density E c that is necessary to satisfy Q = 1. Stable region is below this line and unstable region is above 
this line. The dotted and dashed lines represent the surface densities that is necessary to satisfy Q = 0.5 
and Q = 2, respectively. 

that S c is large enough to satisfy r > 1 in r < r c because of the large Coriolis force and the high 
temperature. Analytic formula for E c is derived from equations (1) and r > 1 part of (11) as 

y _ f _ K A r>a^i6-2b+a -3(7+2a-26)/2 \ 



, KoAB a C^ 2b+a r- i( ' +2a ~ 2b)/2 (r > 1) , (13) 

V256 / 

where A, B, and C are defined in equations (12), (13), and 



7r 2 Gm 

respectively. With the value of k , a, and b in table 1, equation (13) provides us the property 
of E c (see appendix 1 for details). In the fiducial case, in 1.2AU < r < 10AU, sublimation of 
metal dust mainly contributes to the opacity, and S c oc r -171 / 104 . In 10AU < r < 16AU, the 
opacity is dominated by metal dust, and S c oc r" 3 . In 16AU < r < 24AU, the main component 
of opacity is sublimation of ices, and E c oc r~ 7 / 4 . These properties can be seen in figure 2 from 
the fact that the E c line breaks at the place where the main component of opacity changes. In 
r > r c , S c is given by the isothermal condition T = T m i n and Q = 1 as 

= 21 2 ( n -v* " 2 w 1/2 (15) 

ttg s/ viooau; viok; \m q J v ; 

where c Siin i n is the sound speed for T = T min . 

In figure 2, the critical surface density S c changes discontinuously at r = r c . This is 
because the Q value is independent of the surface density there. The reason is explained as 
follows. If the surface density increases, the effect of self-gravity increases to destabilize a disk. 
On the other hand, in the optically thick case, the surface density also has an ability to block 
radiative cooling. Then, the effect of thermal pressure increases to stabilize a disk as the surface 
density increases. If these two effects are balanced, the value of Toomre Q becomes independent 
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of the surface density. This balance is realized when the opacity depends on temperature as 
k oc T 2 in optically thick regime (see appendix 1). The conditions to realize this balance are 
satisfied in the shaded area in figure 1 where ices become the main component of opacity in 
optically thick regime. Thus, the S c discontinuously changes at the critical radius r c where ices 
form. The analytical formula for critical radius r c is derived by equations (1) and r > 1 part of 
(11) with ices opacity as 



In equation (16), r c is an increasing function of a and M s . This is because the viscous heating 
rate is large with large a and large M s . However, dependence on a and M s is weak. Thus, r c 
varies only severalfold even if a or M s is different from typical value by an order of magnitude. 

In figure 2, it is seen that the S c become multivalued function at r ~ 1.5AU and S ~ 
2 x 10 5 g/cm 2 . In this area, molecules are the main component of opacity, and there is the 
inner most radius of the unstable region around this area. We explain this inner most radius 
in appendix 2. 

Note that the critical surface density is derived with the constant viscous parameter 
a. It is considered that this assumption corresponds to the case that the viscosity originates 
from not the self-gravity but the MRI turbulence. This case assumes that the disk is MRI 
active everywhere. Thus, this result is specific to the case with the disk that has no region 
where the MRI is inactive (generally called "dead zones"). In other words, it is assumed that 
there is the significant heating of non-gravitational origin at all radii in the disk. On the other 
hand, the previous studies like Clarke (2009) use the viscous parameter determined by the 
self-regulated gravitoturbulence with the condition Q = 1. In this case, the disk is assumed to 
have the large dead zones without effective viscosity owing to the MRI turbulence. The realistic 
protoplanetary disks must be bracketed by these two extreme cases. Considering such disks is 
beyond the scope of this paper and is left as the future work. 

4. Instability Condition for Protoplanetary Disks 

In this section, the condition for gravitational instability is discussed by comparing the 
surface density of particular disks to S c . We introduce two simple semi-analytic models for 
protoplanetary disks. One is the steady state accretion disk (the steady model). This model 
is realized if angular momentum is sufficiently transported. The other is the disk that has the 
same angular momentum distribution as the cloud core before it collapses. In this model, we 
assume the conservation of angular momentum distribution, and hence this model is hereafter 
called the AMC (Angular Momentum Conservation) model. The AMC model is realized if 
angular momentum is not sufficiently transported. We study the steady model in section 4.1. 
The AMC model is investigated in section 4.2. The interpretation of these two models is 
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Fig. 3. Surface densities with M = 10~ 6 (long-dashed line), 10~ 7 (short- 
dashed line), and 10 -8 M Q (dotted line). Solid line represent E c . We see 
the case of M = 1(T 6 and 10~ 7 M Q have the unstable region at r > 24AU. 



discussed in section 4.3. 

4-1. The Steady State Accretion Disk Model 

First, we consider the model of the steady state accretion disk. The surface density 
profile is calculated in the framework of the standard disk model with a viscosity (Shakura & 
Sunyaev 1973; Pringle 1981). The disk structure are determined by three parameters; viscous 
parameter a, protostar mass M s , and mass accretion rate M. By the equation of continuity 
and angular momentum conservation, we have the relation 



M 
3tt 



(17) 



This equation shows that surface density S is determined by mass accretion rate M. We calcu- 
late the temperature in the same manner as in section 2 using the cooling function approximated 
as equation (8) and opacity represented by equation (10) with table 1, and introducing T min . 
Using equations (3), (4), (5), (11), and (17), we obtain the analytic profile for the surface 
density in this model as 
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where we use a symbols A, B, D, and X are defined in equations (12), (13), 
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Fig. 4. The contours of minimum radius of unstable region rsD.min for various values of 
M s and M. The dashed lines show contours of minimum radius T*sD,min- The area 
where mass flux is less than 7.7 x 10 -8 M Q /yr (left of solid line) has a stable disk. 
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respectively. This steady state accretion disk does not have the inner and the outer radius. In 
other words, the disk is infinitely extended. 

From equation (18) and table 1, we can calculate the Egr> Figure 3 shows the surface 
densities Eg^ for the case with M = 10~ 6 (long-dashed line), 10 -7 (short-dashed line), and 
lO _8 M /yr (dotted line) for the fiducial case defined in section 2. In figure 3, it is seen that 
the lines of Ego break at the radii where the main component of opacity changes just like the 
E c line in figure 2. It is also seen that Esd is large with large M. The critical surface density 
E c derived in section 3 is superimposed in figure 3 by a solid line. The disk is unstable if the 
condition Esd > E c is satisfied. From figure 3, it is found that the disks with M = 10~ 6 and 
lO _7 M /yr are unstable in r > r c = 24AU, where r c is the critical radius defined in equation 
(16). On the other hand, the case with M = lO _8 M /yr is expected to be stable in all radius. 
It is also seen that all the lines including the E c and Esd has the same dependence on radius 
in the outer region (r > 100AU). From equations (15) and (18) for T = T m ; n , we can confirm 
that Esd and E c has the same dependence on radius owing to isothermal state. By using this 
property, we can recast the instability condition Esd > S c as 



M > M crit = 



3acr 
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(20) 



Note that this instability condition (20) is independent of protostar mass M s . It is also found 
that M crit is increasing function of a and T min . As shown in figure 3, the minimum radius of the 
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Fig. 5. Schematic representation of the AMC model. A rigidly rotating cloud collapses to make 
a protostar and a disk. The cloud core is divided into two parts. A protostar is made from 
the part near the rotational axis. A disk is made from the part that docs not become protostar. 

unstable region defined as rsD.min equals to the critical radius r c . The critical radius r c depends 
on the protostar mass as r c oc M s 1//3 , and r c is independent of the mass accretion rate. Thus, in 
the steady model, mass accretion rate determines whether or not the disk becomes unstable, 
and protostar mass affects where the disk becomes unstable. This property described above is 
summarized in figure 4. Figure 4 shows the contours of the minimum radius of the unstable 
region rsD.min in the M — M s plane. It is seen that the stable disk is present in the left area 
where mass accretion rate is low and that the minimum radius is large with large protostar 
mass. It is also seen that the disk becomes unstable at smaller radius than r c for the large mass 
accretion rate M > 5.4 x lO~ 6 M /yr. This property appears where the surface density Ssd is 
large enough to satisfy Ego > E c in r < r c . From figure 3, it is found that Esd is required to be 
larger than 7.3 x 10 2 g/cm 2 in order to become gravitationally unstable in r < r c . 

4.2. The AMC Model 

Next, we consider the disk that has the same angular momentum distribution with the 
cloud core before it collapses. This model corresponds to the disk that has just formed. 

Here, we introduce a simplified model for the formation of a disk and a protostar. Figure 
5 schematically shows the formation process of the disk and the protostar in this model. It 
is assumed that a rotating cloud core collapses to make a protostar and a disk. We divide 
the cloud core into two parts. The shaded part near the rotation axis is assumed to become 
the protostar. The disk is assumed to be made from the other part. The boundary radius 
r cy,in between these two parts is determined so that the mass inside the shaded part equals 
to a given protostar mass M s as a parameter. The cloud core is assumed to rotate rigidly. 
The angular velocity of the cloud core Qq is determined by introducing the rotation parameter 
Po = E'rot/I-E'gravl, where E mt is the rotation energy and E gr!W is the gravitational energy. Note 
that the disk formed in this model has the inner radius r^- m and the outer radius r^out, different 
from the steady model described in section 4.1. It is assumed that the protostar and the disk 
forms instantaneously and that the disk has the Keplerian rotation velocity with ignoring the 
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self-gravity. We also assume the inviscid and axisymmetric formation of the disk. In other 
words, it is assumed that the relation between the mass and the angular momentum is kept 
without redistribution. As the former cloud core, two types of the density distributions are 
considered with the assumption of spherical symmetry. One is the Bonnor-Ebert sphere (Ebert 
1955; Bonnor 1956), and the other is a sphere of uniform density. In this study, we consider 
both cases, and the result for the case with Bonnor-Ebert sphere is mainly described. We 
discuss the difference of these two cases later. 

The Bonnor-Ebert sphere is an isothermal sphere in which the thermal pressure balances 
the self-gravity and the external pressure. This sphere is identified by three quantities, the cloud 
radius Re, the sound speed c Si o, and the central density p c . It is known that a Bonnor-Ebert 
sphere is unstable when the condition 

Re > Rent = 6.46 (21) 

is satisfied. The Bonnor-Ebert sphere whose radius Re equals to -R cr i t is called the critical 
Bonnor-Ebert sphere, and we use this critical Bonnor-Ebert sphere in this model. We set 
c Sj o = 190m/sec, which corresponds to the typical sound speed in a molecular cloud. The 
central density is fixed as p c = 7.45 x 10 _19 g/cm 3 . With these parameters, the critical Bonnor- 
Ebert sphere has the mass M c = 1M and the radius Re = 1.0 x 10 4 AU. Here, we approximate 
the density profile of the Bonnor-Ebert sphere (Tomide 2011, private communication) as 

p(R) = P - m , (22) 

where R is the radius in the spherical coordinate and R c = c S) o/^1.5Gp c . The equation (22) 
approximates the density profile of a Bonnor-Ebert sphere to the accuracy of the error within a 
few percent. This approximation enables us to calculate the surface density of the Bonnor-Ebert 
sphere S B E( r c y ) analytically as 



S BE (r cy ) = 2 / [Re ^' p(R)dz 



\ 1/2 

B- r cy) , 2p c 



D2 2 

(23) 



l + rZ y /R 2 c \{l + R E /R, 
where r cy is the radius in the cylindrical coordinate. 

In this model, the surface density of the disk £AMc(?"d) is determined from the angular 
momentum distribution of the former cloud core as the follows. Once the rotation parameter 
/So is given, the angular momentum distribution of the cloud core is determined. Next, give the 
ratio of disk to protostar mass M&/M s as a parameter, and the specific angular momentum dis- 
tribution of the disk is determined. Suppose that a fluid element in the cloud core at cylindrical 
radius r cy falls onto the disk at cylindrical radius r&. The disk is supported by centrifugal force 
and rotates Keplerian velocity at r&. By the conservation of angular momentum, the relation 
between r cy and is found as 

r cy = (GM s r d /^)V4 . (24 ) 
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r[AU] 

Fig. 6. Surface density with /3 = 0.005 (long-dashed), 0.001 (short-dashed), 0.0002 (dot- 
ted) and Md/M s = 0.25. The critical surface density S c is superimposed by solid line. 



By using the relation (24) with the relation of mass conservation 

r d £ A Mc(Vd)dr d = r cy T, BE (r cy )dr cy , 
we derive the surface density of the disk Eamc as 



J AMC 



(r d ) 



l-r c 2 y /i? 



V2 2 

cy 



(25) 



(26) 



21 + ry^ Vl + i^/flgy r 2 • 
Note that the r cy depends on r d as in the equation (24). From equation (26), we can find that 
Samc depends on radius as rj 2 in the region where r cy > R c . It is also found that Eamc depends 
on radius as r^ 1 ' 5 for r cy < R c . Thus, the enclosed mass in the disk does not diverge. In this 
study, we concentrate on the Keplerian rotating disk without self-gravity. In order to use this 
assumption with validity, the mass of the disk should be less than that of the protostar. If we 
request that the disk rotates Keplerian velocity, the boundary radius r cy in shown in figure 5 is 
always larger than R c . Thus, the Eamc depends on radius as r^ 2 in this paper. By the way, 
realistically, it is considered that a protostar does not have such a large mass just after the 
formation of the protoplanetary disk. In this model, we consider only the outer region of the 
disk and regard the inner region of the disk as a part of the protostar, at least, in the sense of 
gravitational field. The boundary radius r^ in between the inner region and the outer region of 
the disk is determined by the equation (24) and r cyjin . 

Now, we can calculate the surface density of the disk for given M&/M s and (3q. First, 
we look over the property of Eamc f° r the case with fixed M^/M s . Figure 6 shows the surface 
densities Eamc for the cases with O =0.005 (long-dashed), 0.001 (short-dashed), and 0.0002 
(dotted) with M d /M s = 0.25. In the case with large (3 , it is seen that Eamc is small and that 
the disk is present far from the protostar. The solid line in figure 6 represents the critical 
surface density S c derived in section 3. In the case with j3 = 0.0002, the disk is stable because 
the stability condition S c > Samc is satisfied owing to large S c in r < r c = 22AU, where r c is 
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Fig. 7. The unstable condition on Pq-M^/Ms plane for the case with Bonnor-Ebert sphere. The areas la- 
beled as "stable" and "inner" represents a stable disk in the AMC model. The areas labeled as "a few AU" , 
"r c ", and "<30AU" have an unstable disk. "+" symbols denote the cases of figure 6. See text for details. 

the critical radius denned in equation (16) with M s = 0.8M Q . In the case with (3q = 0.001, the 
stability condition S c > Samc is violated in r c = 22AU < r < 54AU, and thus the disk is unstable 
there. This is because S c changes discontinuously at r = r c . In the case with (3q = 0.005, the 
angular momentum is so large that the inner radius r- m = 85AU is larger than the critical radius 
r c = 22AU. The disk can be unstable around the inner radius in r in = 85AU < r < 1.5 x 10 2 AU. 
Note that, in the cases with (3q = 0.001 and 0.005, the outer regions of the disks are stable 
because Eamc rapidly decreases as radius increases. 

Next, we calculate Eamc for various values of M^/M s as well as Pq. Figure 7 summarizes 
the results in (5q-M^/M s plane. The symbols with "+" correspond to the cases in figure 6. The 
results are classified into five groups in the (3 — M&/M s plane. The area labeled as "stable" in 
figure 7 indicates that the disks with these parameters are stable in all radius. The case with 
/?o = 0.0002 in figure 6 belongs to this area. These disks have small outer radius r out < r c owing 
to the small angular momentum. Thus, the condition Eamc < ^ c is satisfied due to the large S c 
in r < r c . In figure 7, the areas labeled as "a few AU" , "r c " , and "<30AU" indicate the unstable 
disks. The area labeled as "a few AU" exists in the small Pq and large M^/M s area. The disks 
with these parameters have small radius and large surface density enough to satisfy Eamc > S c 
in r < r c . Thus, they have a possibility of fragmentation at a few AU from the protostar. The 
area labeled as "r c " exists around middle Pq area in figure 7. This area includes the case with 
Po = 0.001 in figure 6. Although the disks with this area are stable inside the critical radius 
r c due to small £amc> they are unstable in r > r c . This is because S c drastically changes at 
r = r c . It is the formation of ices that enables the disks to become unstable. The area labeled 
as "<30AU" exists at large Pq area. This area includes the case with p = 0.005 in figure 6. In 
this area, the inner radius r in of the disk is larger than the critical radius r c , and the disk is 
unstable around the inner radius. The area labeled as "inner" in figure 7 is difficult to conclude. 
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Fig. 8. The unstable condition on [3o-Ma/M s plane for the case with the uniform sphere. The ar- 
eas labeled as "stable" and "inner" represents a stable disk. The disks with the parameter ar- 
eas labeled as "a few AU", "10-30AU", and "< 30 AU" become unstable at the labeled radius. 

In the present AMC model, the disks with this area are regarded as stable because Eamc is 
very small due to the too large inner radius r- m > 10 2 AU. However, this large r in is artificial 
because the AMC model is too simplified in order to calculate analytically. Realistically, it is 
expected that the disk is present in r < r in . If the disk has the surface density represented as 
in equation (26), the disks with the area labeled as "inner" in figure 7 have a possibility to be 
unstable. 

From above results, we found that the disk is expected to be unstable with large /?o in 
the AMC model. Here, we estimate the critical rotation parameter /3 , c that is necessary to 
become unstable for the disk with small mass (~ 0.01M S ). In this model, the disk becomes 
unstable if the disk extends to the radius r > r c . The outer radius of the disk can be estimated 
by substituting r cy = i? E in equation (24). On the other hand, the r c is represented in equation 
(16). By using the relation 0q oc Qq, we derive the critical rotation parameter /3q,c as 



Note that the r c depends on the viscous parameter a and on the protostar mass M s and that the 
Re and M c depend on the central density p c and sound speed c S) o- In figure 7, it is seen that the 
disk becomes unstable for the flo ^ 0oc even with small mass ratio M d /M s < 0.1. The minimum 
ratio of disk to protostar mass satisfying the unstable condition is Md/M s | m i n ~ 0.0092 with 
0o = 5.0 x 10~ 4 , which belongs to "r c " area. 

We also perform the calculation for the case with the uniform-density sphere as the 
density distribution of the former core. In this case, the core is assumed to rotate rigidly in 
common with the case with the Bonnor-Ebert sphere. This uniform sphere is characterized by 
two parameters, a = E th /\E grstv \ and /3 = i? r ot/|-Egrav|, where a is the thermal parameter and 
E t h is the thermal energy. The thermal parameter is assumed to be ao = 0.86, which is the 




(27) 
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same value as the case with the critical Bonnor-Ebert sphere. The calculation is performed with 
various values of flo and M^/M s . Figure 8 summarizes the results. The results are qualitatively 
same as the case with the Bonnor-Ebert sphere but quantitatively different at some points. 
One is that the critical rotation parameter /?o, c - For the case with the uniform sphere, the value 
Po,c = 1.0 x 10 -3 is larger than that with the Bonnor-Ebert sphere by factor 2.5. This difference 
arises from the fact that the uniform-density sphere has 1.2 times smaller radius Re and almost 
same angular velocity Qq & s compared to those of the critical Bonnor-Ebert sphere with the 
same values of «o and (3q. Since a sphere with a small radius makes a disk that has a small 
outer radius with a fixed angular velocity, the case with the uniform sphere needs larger (3 than 
that with the Bonnor-Ebert sphere in order to satisfy the condition r out > r c . Next, for the 
same radius r in the disk, the mass ratio M&/M s that is necessary for the instability is smaller 
with the uniform sphere than with the Bonnor-Ebert sphere by nearly one order of magnitude. 
In the AMC model, the disk is made from the outer part of the sphere, and the density of the 
outer part with the uniform sphere is larger than that with the Bonnor-Ebert sphere. Thus, if 
both spheres make the disks in same size, the disk made from the uniform sphere has a larger 
surface density than that from the Bonnor-Ebert sphere. Hence, for the same radius r, the disk 
made from the uniform sphere needs smaller M&/M s to become unstable than that from the 
Bonnor-Ebert sphere. Finally, the instability radius around the mid /?o is different. In the case 
with the Bonnor-Ebert sphere, the instability radius is clearly divided into two parts labeled as 
"a few AU" and "r c " . On the other hand, the case with the uniform sphere, the disk becomes 
unstable at the intermediate radius 10AU-30AU. This is due to the large surface density and its 
strong dependence on the radius of the disk from the uniform sphere. All of these differences 
between the disk from the Bonnor-Ebert sphere and from the uniform sphere originate from 
the difference in the relation between mass and specific angular momentum (known as m — j 
relation). Moreover, these differences strongly depend on the assumption of rigid rotation. 

4-3. Interpretation of the Steady Model and the AMC Model 

In previous two subsections, we studied the gravitational stability of protoplanetary 
disks by using two simple semi-analytic models. Here, we interpret these models by considering 
the formation scenario of a protoplanetary disk. 

It is believed that a cloud core with angular momentum collapses to form a protostar 
and a protoplanetary disk. First, the central part of the cloud core begins to contract and 
makes a protostar and a small (typically a few AU) disk (e.g. Bate 1998). Next, the outer 
envelope of the core mainly falls onto the disk, and the mass accretion from the envelope makes 
the disk grow up in mass and size as typically larger than 100 AU (e.g. Adams et al. 1988). 
The time scale for this disk growth is estimated as t grow ~ r j c s .o according to the self-similar 
solution for the collapse of the rotating isothermal cloud (Saigo & Hanawa 1998). This disk 
that has just formed corresponds to the AMC model introduced in section 4.2. Later, this 
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disk experiences the viscous evolution with the transport of angular momentum and eventually 
forgets the initial distribution of the angular momentum. The time scale of viscous evolution is 
estimated as t v i S ~ a~ 1 (H/r)^ 2 Q~ 1 from the azimuthal component of the momentum equation. 
It is expected that a steady state disk, which is discussed in section 4.1, corresponds to the 
state after the viscous evolution. For the case with c Sj o = 190m/sec and H/r = 0.1, the time 
scales for disk growth tg row and for viscous evolution t v - ls are estimated as 

^-"^Mo^liXuf (28) 

and 

^"~ 25yr (i9^7i) '(su)' (29 > 

From equations (28) and (29), it is seen that the disk growth time is much shorter than the 
viscous evolution time for r > 1AU. Thus, it is expected that the viscous redistribution of 
angular momentum can be regarded as negligible during much longer time than the dynamical 
time scale after the formation of the disk. From equation (28), it is found that the viscous 
evolution time is short for a small radius. Hence, it is likely that the steady state accretion 
disk forms at inner region first, and it will spread outward in viscous time scale. In summary, 
we interpret the AMC model as the outer part of the young disk and the steady model as the 
old disk. 

5. Discussion 

5.1. Thermal Stability of Equilibrium State 

We calculated the equilibrium temperature based on the approximated opacity given 
in equation (10) and table 1 (Bell & Lin 1994) in section 2. This opacity depends on the 
temperature, and the property of the opacity varies when the main component of opacity 
changes. Thus, the cooling rate also depends on the temperature similarly. If the temperature 
dependence of the cooling rate is weaker than that of the heating rate, this state is thermally 
unstable (Field 1965). In order to check the stability of the thermal equilibrium state we used, 
we compare the exponent of the temperature in the heating rate to the exponent in the cooling 
rate. In our model, as shown in section 2, the heating rate has the relation T oc T, regardless 
of the optical depth. 

In optically thick case, the temperature dependence of the cooling rate is A cx T 4 ~ b+a ^ 2 . 
Thus, the condition for thermal stability is represented as a — 2b> —6. This condition is satisfied 
except for the case that hydrogen scattering is the main component of the opacity. This case 
is realized where the temperature are larger than the critical temperature T > Tti ~ 3300K for 
p = 10~ 10 g/cm 2 . Within the parameter range we are interested in, the equilibrium temperature 
is safely smaller than Tti, and thus, in optically thick case, the equilibrium state is stable. 
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In optically thin case, the temperature dependence of the cooling rate is A oc T A+h ~ a ' 2 . 
The condition for thermal stability is obtained as a — 2b < 6. From the values in table 1, this 
stability condition is violated when the main component of the opacity is sublimation of ices 
(Ti = 166.8K < T < 202.6K) or sublimation of metal dust (2286.7p 2 / 49 K < T < 2029. 7p^ 81 ). 
According to results in section 2, in the optically thin case, the equilibrium temperature is 
sufficiently smaller than T\ above which sublimation of ices becomes the main component of 
the opacity. Thus, in optically thin case, the equilibrium state satisfies the condition of thermal 
stability. Therefore, the equilibrium state we used in this study is thermally stable. 

5.2. Comparison to Radiation Hydrodynamical Simulations 

Here, we compare the results of our model to the results of some radiation hydrodynami- 
cal (RHD) simulations. First example is the simulation by Boley et al. (2006). They calculated 
the disk evolution about 4000 yr, and the surface density of the disk at the final state was 
approximated as £ ~ 150g/cm 2 x 10~( r / Rjs ) 2 ; where Re = 46.7AU. This disk did not fragment 
in their simulation. Our model also predicts that this disk cannot fragment because this surface 
density satisfies the stable condition £ < S c . Thus, our model is consistent with the result of 
their simulation. 

Next, we make the comparison to the result by Meru & Bate (2010). In their study, the 
evolutions of disks were calculated for the case with different values of the opacities. Their result 
showed that the opacity that is smaller than the interstellar values promotes the fragmentation. 
In our model, the values of £ c and r c represented in equations (13) and (16) are small with the 
small opacity. This means that the disk easily becomes unstable with small opacity. In this 
sense, the prediction by our model is consistent with them. 

Finally, we compare the result of the simulation by Stamatellos et al. (2011). Their 
calculation was performed with nine initial conditions, and they obtained the result that only 
three of nine disks fragmented. From this result, they suggested that the fragmentation is 
unlikely to occur when the disk satisfies the condition that r out < 100AU and that M&/M s <0.36, 
where r out is the outer radius of the disk and M&/M s is the ratio of disk to protostar mass. 
Our result is consistent with their result in the sense that fragmentation is difficult near the 
central star (r < 20AU). However, with their initial surface density profile, our model predicts 
that all of nine disks used in Stamatellos et al. (2011) become gravitationally unstable for the 
case with a < 0.6. These differences between their results and our predictions may be caused 
by the heating due to the gravitoturbulence. In this sense, we should notice that the instability 
condition investigated in this paper is not exactly the same as the fragmentation condition. 

5.3. The Steady Disk Model 

In section 4.1, we derived the critical mass accretion rate M CT \ t that is necessary in order 
that the steady disk becomes unstable. Hartmann et al. (1998) derived the mass accretion rate 
of T-Tauri stars in Taurus and Chamaeleon I, based on the optical observation. The typical 
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value of the accretion rate is as small as M ~ 10 _8 M Q /yr. Since this value is less than the 
critical mass accretion rate M cr - lt ~ 7.7 x 10~ 8 M Q /yr given in equation (20), the disks around a 
T-Tauri stars are expected to be stable. However, in Hartmann et al. (1998), 7 of 40 T-Tauri 
stars in Taurus and 1 of 16 in Chamaeleon I are estimated to have a large mass accretion rate 
enough to become gravitationally unstable. Thus, these statistics indicate that about 18 % of 
the T-Tauri stars have the possibility to have an unstable disk in Taurus, whereas the fraction 
is less than 7 % in Chamaeleon I. By the way, around a protostar, the Keplerian disk has not 
been clearly observed yet. If such a disk is present, its mass accretion rate is estimated as 

M~-L~-±~2x lO' 6 M /yr (_) , (30) 

where Mj is Jeans mass, and tg is free-fall time. Because this mass accretion rate is much 
larger than M cr ; t , the disk around a protostar is expected to be gravitationally unstable. In 
summary, it is expected that most of the disk around a T-Tauri star is stable and that the disk 
around a protostar is likely to fragment. Anyway, the relation between mass accretion rate and 
gravitational instability may be a useful tool to probe the fragmentation in the disk systems 
observationally. 

Clarke (2009) (hereafter C09) also discussed the fragmentation of a steady state accretion 
disk. Assumptions set in C09 are similar to our model, but there are critical differences. C09 
assumes the condition of gravitoturbulence Q = 1, uses spatially- variable a parameter, and 
adopts Gammie criterion as fragmentation criterion. On the other hand, we assume not to 
be in gravitoturbulent state, use constant viscous parameter a, and adopt Toomre criterion. 
These differences produce qualitatively different results. For example, in C09, fragmentation is 
predicted even with very small mass accretion rate (M ~ 10~ 8 M o /yr), whereas, in our model, 
the disk is expected to be stable with such a small mass accretion rate. However, note that both 
studies imply that the formation of ices is an important process. In the regime of ices opacity 
k oc T 2 , the temperature dependence of cooling rate becomes A oc T 2 , which is weaker than 
that of other opacities (sublimation of ices, metal dust, and sublimation of metal dust). This 
property enables the disk to become low temperature when the heating rate becomes small. 
The importance of this fact does not change regardless of fragmentation criteria. 

5.4. The AUG Model 

In the AMC model introduced in section 4.2, the disk is unstable when the condition 
A) > A),c ~ 0.001 is satisfied. By the way, according to the results of radiation hydrodynamical 
calculations, the disk made from the core with the large rotation parameter (3 > /3 ,b ~ 0.01 is 
expected to fragment before the formation of the protostar (Bate 2011). From the observation 
of the line emission of NH 3 and N 2 H + , the rotation parameter /3 of the cloud core is estimated 
to range in 10" 6 < O < 0.07 with the typical value (3 ~ 0.02 (Caselli et al. 2002). Only the 3 of 
20 cores have the rotation parameter /3 < /?o,c ~ 0.001. Thus, only 15 % of the cores are expected 
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to make a stable disk. The cores with middle rotation parameter /3 , c ~ 0.001 < /3q < /?o,b ~ 0.01 
have the possibility to fragment after the formation of a protostar. The 8 of 20 cores is present 
within this parameter range. The number of the core with /3q > /3o,b ~ 0.01 is the 9 of 20 cores. 
The cases with such a large /3q are outside the scope of our AMC model. 

Matsumoto Sz Hanawa (2003), Machida, Inutsuka, Matsumoto (2010), and Tsukamoto & 
Machida (2011) calculated the formation process of the disk from the cloud core using barotropic 
equation of state for various values of the rotation parameter (3q. These studies showed that 
the disks become gravitationally unstable for the cases with large (3q. This feature is consistent 
with the result of our AMC model. Matsumoto & Hanawa (2003) calculated the cloud collapse 
with the rotation parameter 8 x 10 -4 < /3q < 8 x 10~ 2 . They showed that the disk made from the 
cloud with /3 — 8 x 10~ 4 become gravitationally unstable but not fragmenting until the final 
state of their calculation. Our AMC model also predicts that the disk with (3 ~ 8 x 10~ 4 is 
unstable. They also showed that fragmentation occurred before the formation of the protostar 
for the case with /3 ^ 2 x 10~ 3 . These cases are outside the scope of our AMC model. The 
result in Machida, Inutsuka, Matsumoto (2010) is the quantitatively different from our model. 
For example, the value /?o, c ~ 3 x 10~ 5 in their calculation is smaller than that in the AMC 
model by one order of magnitude. However, it is difficult to compare the result of the AMC 
model to that of their calculation quantitatively because, in their calculation, fragmentation 
tends to occur during the phase that the mass of protostar is smaller than that of the disk. 
Tsukamoto & Machida (2011) is more consistent with our AMC model. Despite that they used 
the barotropic equation of state, which is different from us, the value /?o, c = 3 x 10 -3 in their 
calculation is consistent with our AMC model within the error of factor 3. 

5.5. The effects of ignored processes 

In this study, we ignored some physical processes that are expected to affect the grav- 
itational instability in the disk. Here, we comment on the effects of processes ignored in this 
paper. 

First, the self-gravity will affect the disk properties through the effect on the transport of 
angular momentum and energy dissipation. In the form of gravitational torque due to the non- 
axisymmetric mode and the gravitoturbulence, the self-gravity affects the viscous parameter 
a. It is discussed that the a becomes large if the Q value approaches unity (Kratter et al. 
2008). However, this a as a function of disk properties has not been clarified yet. In this study, 
in order to avoid this uncertainty, we used constant a parameter and derived the condition 
for gravitational instability as a function of a. The work considering a from the self-gravity 
remains as a future work. The self-gravity will also affect the scale height H. If considering 
the self-gravity in the vertical direction, the scale height H sg with self-gravity is represented as 
a function of Q value. This H sg for the case with Q = 1 is approximately half of the H without 
self-gravity represented in equation (5). From this effect, the values of S c and r c are modified, 



20 



but the amount of this modification is at most a few tens of percent. Thus, it is considered that 
our result does not qualitatively change even if considering the compression by self-gravity in 
the vertical direction. 

Next, the effect of irradiation from the protostar will affect the disk properties. At the 
region far from the protostar, the heating by irradiation from the protostar becomes greater 
than the viscous heating because the viscous heating strongly decreases as the radius increases. 
Thus, at the far from the protostar, there is a possibility that the disk is stabilized by the 
irradiation heating. The heating rate by irradiation is affected by the disk flaring, and this 
flaring is determined by the temperature distribution of the disk. It is difficult to make the 
accurate treatment of irradiation heating in our model because we did not consider the radial 
distribution of the disk when deriving the instability condition in section 3. However, the effect 
of irradiation heating is weak in disks with large surface density. Thus, it is considered that 
our results for optically thick case will not critically change. This irradiation heating should be 
taken into account when investigating the properties of the particular disk. 

Finally, we comment on the effect of magnetic field. The magnetic field has much 
influence on the angular momentum transport during all phases from the protostar formation 
to the disk evolution. In the phase of the protostar formation, it is considered that the magnetic 
field triggers the outflow and that it transports the mass and the angular momentum (Shu et 
al. 1994; Tomisaka 2002; Machida et al. 2008). Thus, it is expected that the assumption used 
in the AMC model is violated by the effect of this outflow. As a future work, it remains to 
construct the model including the effect of the outflow. In the phase of the disk evolution, the 
turbulence arising from the magnet rotational instability (MRI) is considered to become the 
main source of viscosity. In the region where ionization degree is low (called as dead zone), the 
MRI turbulence is inactive (Sano et al. 2000). Thus, in the dead zone, the viscous parameter 
a is expected to be much smaller than that in the MRI active region. Although we used a 
constant a model in this study, as the future work, it will be important to consider the realistic 
viscous parameter based on the physical processes with the relations among other physical 
quantities. 

6. Conclusions 

In this study, the gravitational stability of protoplanetary disks is studied. The temper- 
ature in a protoplanetary disk was analytically calculated with considering thermal effects such 
as radiative cooling, viscous heating, and ambient heating source in a molecular cloud. It is 
found that the temperature becomes as low as T <20K in optically thin regime. We derived a 
critical surface density S c , which is necessary for the disk to become unstable, as a function of 
radius. By comparing a given surface density to E c , the possibility and the position of the grav- 
itational instability can be predicted. The formation of ices is important for the gravitational 
instability, and we found that the S c changes discontinuously at the critical radius r c ~20AU, 
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where ices form. 

Two semi-analytic models of protoplanetary disks are used with above critical surface 
density S c in order to discuss the possibility and the position of gravitational instability. In 
the model of a steady state accretion disk, which corresponds to the evolved disk, it is found 
that the disk becomes unstable when the mass accretion rate is greater than the critical mass 
accretion rate M cr i t in equation (20). From this result, it is expected that most of the disks 
around T-Tauri stars are gravitationally stable and that the disks around protostar are unstable 
in r > r c . In the model of the disk with the same angular momentum distribution as the former 
cloud core, which corresponds to the outer region of the young disk, the disk is expected to 
become unstable in r > r c for the case with large rotation parameter (/3q ~ 10~ 3 ). From the 
results of these two models, we conclude that the fragmentation near the central star (r < r c ) 
is expected to be rare to occur as long as the parameter range indicated by observations is 
considered. 

We thank Fumio Takahara, Yutaka Fujita, and Hideyuki Tagoshi for fruitful discussion 
and continuous encouragement. We are also grateful to Taishi Nakamoto and Kei Tanaka who 
provided helpful comments and suggestions. 

Appendix 1. The expression of Temperature, Q value, and Critical Surface 
Density 

Here we summarize the value and dependence of the temperature, Toomre Q value, and 
the critical surface density S c . 

A. 1.1. Equilibrium Temperature 

The equilibrium temperature given in equation (11) is explicitly written below. The 
temerature range for each formula is given in table 1 
A. 1.1.1. Optically Thick Case 

(a)ices 

t=12K (toau) (lovcmO im) (A1) 
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(c) metal dust 

(d) sublimation of metal dust 
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A. 1.1. 2. Optically Thin Case 

As described in section 2, the temperature become very low in optically thin case. Thus, 
we consider only the case with ices opacity. 

(a)ices 

™(lAu) (lit) (oiii) (A5) 

4J.j8. T/ie Q Fake 

The Q value defined in equation (1) is represented below. 
A. 1.2.1. Isothermal Region 

A. 1.2.2. Optically Thick Case 

In optically thick case, Q value is represented as 

(27 \ 6+a-26 3 7+2a-2b 26-4 

— n AB a \ r -2^+^2irS^+^F, (A7) 

where A, B, and C is defined in equations (12), (13), and (14), respectively. Using the value 
of table 1, we can see the dependence of Q value. From equation (A7), it is found that the 
Q value is independent of surface density for the case with b = 2, which is realized when ices 
dominate the opacity. 
A. 1.2.3. Optically Thin Case 

In optically thin case, Q value is represented as 

/ 27 A \ 6-a+2b 3 2a-2h-7 -6-26 
Q = C[ r 2 6-a+2b £S^+2S, (A8) 

where A, B, and C is defined in equations (12), (13), and (14), respectively. 
A. 1.3. The Critical Surface Density S c and the Critical Radius r c 

(a) isothermal region 

Cs,min^ 91 , 2 fr_y'^ fT min \^ 2 (M s \ 1/2 

-o-^ = 21 g/cm j j ^— j (A9) 

This is available for r c < r, where r c is the critical radius represented as equation (16). 

(b) ices 

In this regime, unique E c is not present but r c is present. 
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(c) sublimation of ices 



V10AU/ V 1M 0/ ^°- 01 ' 
This is available in t\ < r < r c , where T\ is the radius where the main component of opacity 
changes between sublimation of ices and ices. This T\ is represented as 

M s \ 1/3 / a \ 2 / 9 



n = 16AU 
(d) metal dust 



lMfT 



U) ^ 



Sc = 6.2 x 10 W " 3 f ^ V (—) (M3) 

This is available in r 2 < r < r\, where r 2 is the radius where the main component of opacity 
changes between metal dust and sublimation of ices. This r 2 is represented as 

/ M, \ 48/141 / a \ 98 / 423 
r 2 = 10AU (-^-) (A14) 



s im ; vo.oi 

(e) sublimation of metal dust 

, , „ / r \ ~ 171 /i04 / Ms \ 7/13 / a \ x / 52 , . . 

S c = 6.3 x 10 3 g/cm 2 — — — - A15 

This is available in r e d gc < r < r 2 , where r e d ge is the inner most radius for the gravitational 
instability represented as 

( A'L \ 1/3 / a \ 54 / 353 
^.OAU^j (_) ( A16) 

We explain this r e d ge in Apendix 2 

Appendix 2. The inner most raius for the gravitational instability 

In the regime where molecules mainly contribute to the opacity, which is realized at 
1600 K, the temperature dependence of cooling rate becomes weak as A ocT 4 / 3 . Then, the 
temperature becomes too high with large surface density, and the stabilizing effect of thermal 
pressure becomes larger than the destabilizing effect of surface density increasing. It means 
that increasing surface density stabilizes the disk. Thus, the line of S c become multi-valued 
function of radius. In figure 2, we see this feature for E ~ 2 x 10 5 g/cm 2 and r ~ 1AU. The 
upper line of S c slopes upwards when going from left to right. This is due to the lower heating 
rate at larger r. Because the S c line slopes downwards where the main component of opacity is 
sublimation of dust grains, the unstable region has the inner edge at the point where the main 
component of opacity changes between molcules and sublimation of dust grains. In figure 2, we 
see that this inner edge locates at r ~ 1.2AU and S~2x 10 5 g/cm 2 . We name this radius r ec j ge . 
For r < r e d ge , the disk cannot become gravitationally unstable no matter how large the surface 
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density becomes. By using the equation (11), (13), and the condition for opacity-change from 
dust sublimation to molecules, r e d ge is represented as 



This r cdgc weekly depends on a and M s just like a case of r c . Thus, even if a or M s is much 
different from typical value, the value of r e d ge varies only severalfold. 
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